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ABSTRACT 


A stable numerical procedure for solving the coupled, 
nonlinear equations of electron and ion conservation and 
electrical potential (Poisson's equation), is described. A 
two-dimenSional model with periodic active sites on a flat 
plate is utilized to obtain both qualitative and quantita- 
tive results which clearly illustrate the self-generating 
sheath and ambipolar regions adjacent to a non-emitting 
electrode. The active sites are characterized as voltage 
sources and by electron densities depressed from both the 
non-active wall and the free-stream values. Various cases 
of species density at the nonactive wall are evalnated. 
Recombination/ionization plays an important role in estab- 
lishing the boundary layer behavior in the ambipolar region 
and in the dimensionality of the problem formulation. 
Application is for Nitrogen gas at one amagat in a pulsed 
discharge of about 50 micro-seconds in duration. 

Results clearly demonstrate the boundary layer nattre of 
the species density at the electrode. A sheath length of 38 
microns about a 35-volt active site and 5-20 microns along 
the non-active portion of the electrode is established. 


Joule heating is determined to be not important. 
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I. INTRODUCTION 


A. BACKGROUND 

Losses at the electrode boundary regions play an impor- 
tant role in high energy density molecular gas flow lasers. 
Of particular interest to this study is the Carbon-Dioxide 
gGlow-discharge plasma, of which Nitrogen is the principal 
energy bearing constituent. 

The main energy loss mechanisms in a collision dominated 
plasma are found in the sheath and ambipolar regions. The 
sheath voltage loss occurs directly as a result of Debye 
shielding, which creates a space charge layer abutting an 
electrode or wall. The ambipolar losses are associated with 
ambipolar diffusion, relating to the charged species density 
being less than the free-stream or equilibrium density. 

This loss exists even in the absence of current flow. 

An initial premise of probe theories that takes account 
of collisions between charged particles, atoms and molec- 
ules, is that the plasma is guasi-neutral up to a certain 
distance from the probe and that the electric fields are 
localized within a narrow sheath region adjacent to the sur- 


face of the probe. The problem of making the transition 
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from the quaSi-neutral solution to the sheath solution was 
usually resolved through matching techniques. The descrip- 
tion of a collision dominated sheath together with the qua- 
Si-neutral ambipolar region presents an imposing set of cou- 
pled nonlinear partial differential equations, for which 
there exists no known exact solution. This study reports a 
humerical method that describes the sheath and ambipolar 
regions in a self-consistent fashion. 

In the hydrodynamic approximation (i.e., coilision domi- 
nated), the dimensions of the space charge region "dA," (the 
sheath length) is larger than the mean free path "¢" of the 
species particles. Additionally the sheath characteristic 
length is very auch smaller than the characteristic length 
of the plasma boundary layer. This puts our study in the 
collisionally thin sheath region, {Ref. 1]. 

However, the plasma is not considered to be in thermal 
equilibrium. In fact for plasmas in thermal nonequilibriun, 
the electron temperature is governed by the electron energy 
equation, which includes the effects of the electron energy- 
conduction and electric field as well as collisional energy 
exchange between electrons and heavy-gas particles. The 
consideration of negligible energy exchange during colli- 


sions (thermally frozen) does not generally result ina 
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constant electron temperature (Tg); particularly in the 
sheath region, even if other properties such as the ion 
temperature (T,) are constant. Therefore the electron 
energy equation is necessary even in the case of the ther- 
mally frozen situation. 

An ionizer/sustainer discharge pumped laser requires 
that relatively high E/n values be distributed as uniformly 
as possible throughout the bulk of the gas medium. The 
electrode regions however represent a substantially higher 
and more inhomogeneous E/n which may impose an instability 
on the remainder of the discharge. In the region of the 
electrode, the E/n is significantly higher than the undis- 
turbed plasma. In this region tne alevated electron temper- 
ature, which is a direct function of the sustainer field, 
Will increase ionization, which could lead to arcing and 
breakdown. 

In the so-called vacuum arc [Ref. 2], anode spots appear 
under certain operating conditions. These anode spots are 
of practical importance, having been also observed in high 
density plasmas [ Ref. 3]. They are typically micron or 
sub-micron sized protrusions on the electrode surface, which 
act to "focus' the electric field and current stream lines. 


Primary consideration here is with the anode. 


WZ 








Be. REVIcW of PREVIOUS WORK 
We whiysics 

Most work on Sheath phenomena is embodied in “probe" 
theory investigations, where the affects of the alectrode 
(probe) locally disturb the guiescent plasma. Such work is 
relevant to this study since the anode is essentially a 
heavily biased probe in contact with the plasma, which is 
quiescent within the yeas region. 

Analyses of electrostatic probes in collision doni- 
nated plasmas were originally done by assuming that a quasi- 
neutral diffusion controlled reqion extended to within one 
mean free path of the probe where it was matched to the edge 
of the free fall-Sheath. No provision was made for a tran- 
Sition region between the two regimes [Ref. 4]. 

Shultz and Brown {Ref. 5] studied the transition 
from a collisionless to collision dominated sheath through 
charged species collection with a spherical probe. Postu- 
lating different physical models for cases ranging from a 
collisionless sheath to a sheath with many collisions, 


semi-empirical formulas were obtained in each instance. 





Several other authors attenpted matching the solution in the 
collision dominated region directiy to that of the free-fall 


Sheath uSing variational principles. 
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The previously mentioned techniques can be criticized, 
because in each of these methods something was forced to 
fit. In other words, at least one derivative was disconti- 
nuous at some artificially imposed boundary, thereby raising 
doubt as to the validity of the results. 

The first systematic analysis of probes in collision 
dominated plasmas were carried out by Cohen [Ref. 6] and Su 
and Lam {ref. 7], who assumed the continuum equations valid 
throughout the plasma, including to the probe surface. Asa 
consequence of these analyses, both the sheath and the qua- 
Sineutral regions appeared naturally from the diffusion 
eguations in the limit as the ratio of probe size to Debye 
length approached infinity. The continuum model removed the 
concern about forcing sets of equations to match at the edge 
of an arbitrarily assumed sheath. These authors restricted 
the problem to an isothermal plasma, but as will be shown 
here, the electron temperature may vary significantly in the 
neighborhood of the probe for most discharges of interest. 

Radbill (Ref. 8} yuasilinearized the continuum equa- 
tions, and solved with numerical integration, obtaining 
results that agreed with those of Cohen where comparison 


could be aade. 





In the previous analys2s, not only was the isother- 
mal plasma assumed, but the plasma was so slightly ionized 
that only charge-neutral interactions needed to be taken 
into account. However, for highly nonequilibrium plasmas, 
charge-charge collisions play an important role {[Refs. 9 and 
10). The effects of collisions between charged particles 
manifest themselves as volume ionization, volume rscombina- 
tion and the dependence of the transport coefficients on 
Charged particle density and temperature and, therefore on 
position. Barad [Ref. 11] cites a Specific example, 


"If the electrons have a tenperature of 2500 °K, then the 
lonization fraction egual to about 107-5 will result in 
charge-charge collisions. occurring, one-tenth as frequently 
as charge-neutral collisions. [t is shown... that for 
such a degree of charge-charge interaction, significant 
effects are felt." 


Chung, Talbot and Touryan (Ref. 1], in an extensive 
review of probe studies, state that no general solution is 
available for determining charge density and species temper- 
ature for probes small in comparison to the hydrodynamic 
boundary layer thickness. 

Dolson {[Ref. 12] investigated the nature and extent 
or the voltage drops in the vicinity of magnetohydrodynanic 
(MHD) non-emitting electrodes, in particular the losses 
attributable to the sheath. Under his conditions the nonex- 


istence of a one-dimensional sheath solution 1S shown, and a 
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computer model with two-dimensional, periodic active sites 
representing a flat plate electrode is developed. With this 
model the effects of a magnetic field and Joule heating are 
studied and the results are compared with experimental 
observations reported by Argryoupolos (1973), and Sonju and 
meno (1978), [ Refs. 13, 14]. 

An extension of Dolson's work is reported in 1980 
(Ref. 15], for application to ionizer/sustainer discharges. 
The role of electron pressure is introduced along with net 
ionization and recombination, which oe necessary in 
establishing truely boundary-layer behavior for the charged 
species profiles. Additionally, a modified Newton-Raphson 
computational procedure, which enables much faster solution 
convergence is formulated. That study demonstrates solu- 
tions with non-active wall densities depressed from the free 
stream density, a result Dolson did not accomplish. Results 
were for 10 volts of electrode potential fall, an increase 
from 5.0 volts reported by Dolson. 

2. Computational Methods 

The fundamental difficulty in the development of 
computational models for the prediction of plasma flow phe- 
nomena lies in the coupling of equations governing both the 


gasdynamic flow and the electromagnetic (EM) fields. Since 


a 
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the flow fields can affect the EM fialds (and vice-versa), 
the nonlinear partial differential equations must be solved 
Simultaneously and the formulation of the problem becomes 
complicated. In the present study the flow properties are 
decoupled from the the charged species conservation equation 
and Poisson's Equation. 

The approaches used in the solution of plasma flow 
phenomena are quickly reviewed below so as to point out the 
background of the the approach in this study. This problen 
is typical to those dealing with plasma-boundary Layer phe- 
hogwena. Approaches to solving the general plasma flow prob- 
lems have consisted of the following; 

(a) One dimensional "influence coefficients." 

(b) One-dimensional closed form analytic solutions to 
ordinary differential equations. 

(c) Similarity transformation and closed form solutions 
to a modified differential equation in one- or two- 
dimensions. 

(d) One-, two-, or three-dimensional numerical solution 
to the uncoupled steady plasma flow phenomena 
equations. 

(ce) One-, two-, or three-dimensional numerical solution 


to the coupled steady plasma flow phenomena equations 
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using either the magnetic induction equation or the 
elliptic equation governing the electric fields. 

(tC) One-, two-, or three-dimensional transient numerical 
solution to the coupled plasma flow phenomena 
equations using either the magnetic induction or the 
elliptic equation. 

Methods (a) through (c) are exact but are restricted 
to a narrow range of problems. Approaches (d) and (e) pro- 
vide rapid convergence to a steady state solution but do not 
consider the transient behavior between the gasdynamic and 
electromagnetic fields, or often the strong coupling between 
the two sets of fields. Method (f£) is most desirable since 
it treats a wide range of problems, but suffers from large 
computational times [ Ref. 16]. The approach of this study 


1s characterized as a method (e) approach. 


Co GB@SCTIVES OF 'THE PRESENT STUDY 

Fhe objective of this study is to demonstrate a numeri-~ 
Cal procedure which can model the physical behavior of a 
conducting fluid, as would be generated in the after-glow of 
an electron beam discharge laser, under the influence of 
electric forces. Calculations are performed in two dimen- 
Sions, with a flat plate non-emitting electrode (anode) 


which has periodic active sites. These active sites 
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facilitate current constriction required to satisfy charge 
continuity, Ohm's law, and Poisson's equation. The physical 
model is idealized in Fig. 1. 

The controlling equations are solved for the potential, 
the charged species density distributions, electron tempera- 
ture and current in the vicinity of the electrode. The 
resulting computer programs generate sheath and ambipolar 
regions in a self-consistent fashion, using the same set of 
equations throughout the field. The size of the sheath and 
the voltage drop attributable to its shielding effects are 
investigated. 

This work does not consider those phenomena that lead to 
the arc discharge, which results from thermal instabilities. 
The constriction at the anode is purely that required to 


Satmscy continuity. 
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A. GOVERNING EQUATIONS 

The governing equations that describe the electrode 
region portion of the plasma, consist of Poisson's equation, 
species conservation and electron energy equation, in addi- 
tion to the overall continuity, momentum and energy equa- 
tions of the plasma. These overall plasma equations are 
essentially those of the background gas since the plasma is 
weakly lonized (less than 0.01 4). 


We have [Ref. 1]; 


Poisson equation... 


2 
Boe a : 
V ¢ an £ (ne Zni) (1) 
electron conservation... 


U-Vne - V-| De (EVNe - £- Ke 7) | = Ne (2) 


ion conservation... 
Tr -7-(p (Mon « LK = hi (3) 
UVA:- V | 0; (5:90 + <n vo)| 
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In gas discharges, the electron energy distribution is 
known to be non-Maxwellian. In order to allow for non-Max- 
wellian distributions and still retain use of a temperature 
concept, an empirical fit to the calculated relation of 
electron temperature as a function of E/n is used in the 
representation of the electron temperature (sae Sec. III.A). 
The concept of a "two-temperature" plasma is preserved but 
no restrictions are made on the energy distribution function 
of the electrons. In addition to the assumption of a func- 
tional dependence or the electron temperature upon E/n only, 
there is no representation of an external magnetic field 
and, in particular the internal or self-generated magnetic 
fields are not considered. The overall gas temperature is 
assumed constant and in equilibrium with the plasma ions. 
The bulk energy equation is assumed to have little bearing 
on this study, since there is negligibie gas heating as dis- 
cussed in Section III.A. This effectively uncouples overall 
energy consideration from the problem. 

Tre electrodes serve the purpose of sustaining the E/n 
to pump the laser medium. No account is made of an external 
E-beam or other source of primary ionization. Simply, the 


medium may be considered to exist in an afterglow state. 


22 





The production term (i.e., the right hand side) in Eqns. 
(2) and (3) has been described utilizing measured and theor- 
etical results of the rate coefficients for Nitrogen ioniza- 
tion and two-, and three-body recombination. The production 


term 1s expressed as; 


7 7} De Kn,n; - %&ne2nz- Seles (4) 


where "yy." and "@" are respectively the ionization and 


recombination rate coefficients. 


3. BOUNDARY CONDITIONS 

In the customary poundary-value problem in fluid mechan- 
ics, a two-dimensional flat plate is assumed. A similar 
two-dimensional Cartesian description is used in the present 
work, which has been found to be the mininum suitable 
description for the problem posed, {Refs. 3,12,15]. No 
One-dimensional solution exists for the frozen flow of 
charges through a collisional sheath, since a Cartesian geo- 
metry in one dimension is devoid of the necessary geometric 
decrease of current density away from the electrode, as 
there would be in a cylindrical or spherical geometry. 

Convection is considered negligible, as can be demons- 


trated by comparing fluid velocity with the drift velocity. 
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Electron and ion densities are determined at the plasma 
edge, knowing the state of the medium. A current balance 
can be used to establish the species densities resident at 
the active portion of the anode. An electrode voltage drop 
1s determined from the properties of the medium; e.g., the 
effective ionization energies of the gas in particular. 
Since current constricts, there must be regions of the flat 
plate that receive no current. The inactive portion of the 
flat plate merits attention, especially in the regions adja- 
cent to the active anode site. Species concentration along 
th2 inactive portion of the wall can vary, depending on the 
reactivity (catalytic properties) of the wall. ‘There is a 
Dalance that can be established between the voltage varia- 
tion normal to the wall, and the species concentration at 
the wall together with the species variation normal to the 
wall. A zero concentration of electrons and/or ions at the 
wall is an acceptable boundary-condition; which is to say 
that the wall acts as a perfect sink for the charged spe- 
cies. A detailed kinetic-theory analysis shows that there 
is in fact a lower bound to the electron and ion concentra- 
tions at the wall {Ref. 17]. In any case, there is no cur- 
rent into or out of the inactive wall, (see Pig. 1). This 


ls discussed later in Sect. IV. 
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The remaining portion of the boundary conditions refer 


to the iree stream medium, and are easily established. 


C. PARAMETERS OF THE DISCHARGE MEDIUM 
The application of this study is for molecular Nitrogen 
at one amagat, where the gas temperature is 273 °K. Other 
pertinent parameters follow; 
Meeestic Field: 1.5 to 15 KV/em. 
Zlectron Density: 10+!7 +o 1019 g-3, 
Electrode Voltage Drop: 35 volts. 
The electron energy relation is, as mentioned, a curve fit 
from 2xperimental results [Ref. 18], carried out by RW. 
Compton and D.J. Sutton (1952). Since cross-section data 
for electrons in Nitrogen discharges are available [{ Ref. 
19], the 2lectron temperature and diffusion coefficient are 
readily determined, and representable as a function of E/n. 
The ion temperature is assumed to follow that of the 
neutral gas. An electrode potential fall of about 35 volts 
may be assumed since it is the effective ionization energy 
for Nitrogen gas and, therefore, the anode fall [Ref. 20]. 
The cathode is more complicated to describe because of 
the requirement of surface electron emission. The cathode 
fall in cold cathodes is of the order of 250 volts, and is 


considerably less in thermionic cathodes. A proper 
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III. PROBLEM SOLUTION 


A. PROBLEM SIMPLIFICATION (PHYSICS) 
There are a number of simplifications that are incorpo- 
rated in the problem description. Namely; 
e Steady state conditions (within the pulse time). 
e No magnetic fields. 
e Negligible convection. 
e Negligible Joule heating, 7, = T,= constant. 
However, specific account is made of; 
e Ionization /recombina tion. 
e Electron pressure. 
These are discussed in part in this section. 
1. Convection 
The ambipolar region is the transition region fron 
the sheath to the undisturbed plasma and can perhaps span 
the boundary layer. This fact reguires some consideration 
of the ee of convective effects. 
From the ion species conservation equation (Eq. 3), 
for example, we know that in the ambipolar region the last 


term of the left hand side is small compared to electric 


conduction. Now we need only to compare the convection with 
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the conduction, or simply the fluid velocity with the drift 
velocity. 
v= YR + 0.1 to 1.0 n/sec 
where the Reynolds number (R,) is 10* to 106 and the free 
stream velocity is 10 to 100 m/s. Now the drift velocity 
may be approximated by 
Vy = HES = 10 to 100 o/s 
for E = 105 to 7106 V/ym and the mobility “ = 10-* 
m@/s/volt. Clearly, in the presence of a relatively strong 
interelectrode field, the contribution of ion convection due 
to a cross flow should be negligible. The geometrical 
Orientation of the discharge with respect to the primary 
slow can thus become an important consideration. 
2. Density Boundary Layer Formation Time 

The charged-species density profiles change ina 
more complicated manner than either the voltage potential or 
the electric field. The magnitude of the ion and electron 
densities may change appreciably within the Sheath; more- 
over, fractional analysis is risky because inflections are 
present in the profiles. A discussion pertaining to the 
Stability of the ambipolar region and its boundary layer 
nature is presented in detail by Biblarz, et al., [{Ref. 15]. 


However, Appendix A shows that the characteristic length of 
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the ambipolar region is of the order of the characteristic 
Boemen sength. betting L=5x10-% meters and the drift veloc- 
ity as before, vy = 10 to 100 m/s, a characteristic time for 
the formation of the density boudary layer would be 0.5 to 
5.0 microseconds. It is reasonable that the ambipolar dif- 
fusion to the wall can establish itself well within the 50 
microsecond pulse period of the electrodes. 
3. Significance Of Joule Heating 

For the present problem, joule heating can be shown 
not to be Significant. [In the vicinity of the active site, 
the average energy deposition rate per unit volume is shown 
to be approximately 4x108 watts/m3. The temperature of the 
Nitrogen gas at one amagat would be raised approximately 1°K 
in 2.6x10-6® sec. Allowing a 50x10-® sec pulse period would 
raise the temperature of that region by approximately 209K, 
Gr 7'%- 

Assume that the heat generated at the anode spot is 
Carried away by thermal diffussion (conduction). Then for 
Nitrogen gas (k=0.0267 W/m-°K, fp =1.138 kg/m3, and cp, =1043 
J/kg-°K) the thermal diffusivity is; 

= gos = 2.25x10-3 m@/sec 

A typical diffussion time would be (for, Ay =4.4x10-5 a); 


2 
tye ds/l= 0.86 millisec 


Zo 


ee: eae eae 





The "flush time" is the time required for clearing 
the lasing cavity of residue gases in preparation for the 
next pulse, and can be estimated as the gas flow rate 
divided by the flow length. Assuming v = 100 m/sec and L = 
10 cm, the flush time would be: 

FT = 0.1/190 = 1 millisec 
Other factors considered, the inter-pulse period (time het- 
ween pulses) would be larger, but even in a one millisecond 
period, there is just enough time to cool the anode spot. 

The analysis so far has assumed a flow velocity at 
the surface to be the same as the bulk flow in the laser 
cavity. However a boundary layer exists-- probably turbu- 
lent, with a laminar sub-layer adjacent to the electrode 
Poms 

i> = xV,/) = (0.1) (100) /(15.53x10-%) = 6.44x165 

The velocity boundary layer thickness would be 2.55 milline- 
ters. Assuming a 1/7-th power profile, the velocity ata 
Sheath length from the electroije would be: 

v= Ye (A/S) = 100 (4.4x 1075/2.55x10-3) ¥7 = 55 m/sec 
A characteristic local flush time would be approximately 2 
milliseconds, and much nigher for a laminar velocity pro- 
file. These flush times can be quite large as the require- 


ment to cool the Joule heated gas near the anode spot is 


30 





increased. The design flush time based on other factors may 
infact be insufficient to clear enough of the Joule heated 
gas from the vicinity of the anode spot, enabling a build up 
of heated gas, which could lead to local enhanced ionization 
Oz the gas. 
4. Electron Temperature 
The electron temperature (I,) is assumed to be a 

linear function Of E/n as presented in Fig. 2, for the 
region of interest, namely E/n = 0.5 to 5.0x10-29 V-m@. For 
the range of E/n considered, it is adequate to take the 
electron temperature as; 

Te = [Aelog(E/n) + BJjeT, (5) 
where; 


A eet 


B 3325 
E/n is normalized to 10*2°9 V-n2 
The diffusion coefficients were determined from 
Brown (Ref. 18]. The drift velocities were obtained from 


plotted results, and since (=v /E, and De = MykTs /e, (Ein- - 


stein's relation), the diffusion coefficients are; 


De 0.0665 m@/sec 


D. 


4 


8.0x107-S me/sec 
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Be. @ORKING EQUATIONS 


Since convection is negligible, and assuming that the 
sustainer operates in the afterglow of an electron beam, the 


governing equations, namely Eqs. 1,2 and 3, can be written 


as follows; 


“~ 
A 
V2d = c, (Re- 1g) (6) 
A PA a” A Af ~ 
-Vef Vn, + nye /Q - no V4/6 | =Co Ne (7) 
“ AA AAA “A 
-Yo( Vn; + On; VG] = Can (8) 
where; 

A A 

Ng= Dy /Negp P= ef/kTe, 

OE ©. Te, P= Tey/Te 

he e¢/kTens Cy= AZNey/D, Cz3= C2D,/DE 


The equations may now be applied to flat plate formula- 
tion as depicted in Fig. 3. It is convenient to drop the 
"a" for simplicity, but it is to be understood that all 
variables have been suitably non-dimensionalized. 

de have Poisson's equation; 


od a qo = ¢,[ne-ni] 


Fx oe (9) 
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electron conservation equation; 


Cy Ne(Ne-ni) + one 26 ane 9d] _ Ne(ddde sf 22 | 
= Pa "Sy oy 8 | ax ox Dy oy 


“| 3 Ine . dO Ire a | an ne. 2 st | ros 


Ox OX oy Oy |" 
+ Mel (@ (26+ Bey | - ne| 28 x +38, = C, One 


Hon conser vation; 


CAN: (te-Ni) “6 [3 a + an 34} - -[ Sas + dae] = Cen 


C. NUMERICAL METHODOLOGY 

The nonlinear terms creat problems in the computer anal- 
ysis of the electron and ion species conservation equations, 
(Eqs. 10,11). The Jacobi method includes all nonlinear 
terms on the "right hand side"; i.e., external to the coef- 
ficient matrix. Convergence to a solution 1S possible if 
these non-linear terms change slowly over each iteration. 


It was found by Dolson [Ref. 12] that the Jacobi method was 
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in fact unstable for the present set of equations and 
conditions. As a consequence of the failure of the Jacobi 
method, a quasi-Jacobi method was implemented, in which an 


estimate for each of the solution variables (F, A n-) was 


ee 
computed. When the product of two variables is encountered, 
one variable is treated as a constant coefficient for each 
iteration. This means that the non-linear terms are 
retained in the coefficient matrix. The "constant" coeffi- 
Clients are updated after every iteration, thus altering the 
coefficient matrix. The conventional Jacobi method was 
found to converge only for low values of voltage potential 
at the electrode, whereas the quasSi-Jacobi procedure pro- 
vided converged solutions for approximately 5 volts. A dis- 
cussion of these methods is provided in Appendix ¢ as back- 
ground for the development of the modified Newton-Raphson 
(MNR) method, which is used in the solution to Fqs. 1 
through 4. 

The nonlinear coupled partial differential equations 
(PDE), Eqs. 9,10 and 11 are linearized to provide a set of 
linear algebraic equations, which must hold throughout the 
domain described as follows. The domain as illustrated in 
Pig. 3, is divided up into equal sided sub-areas whose ver- 


tices are the computational nodes at which the linearized 
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algebraic equations must hold. The wall is represented by 
seven (7) nodes, with the leftmost node being the location 
of the active site. The "y" direction extends from this 
active site along the wall. The "x" dimension perpendicular 
to the wall is represented by 25 nodes starting from each of 
the wall nodes. The solution to the PDE's is represented by 
the vector 2Z(x,y), or on the eee ional mesh aS Z(i,j)- 
The boundary conditions have been discussed praviously 
and are presented in Pig. 3. The boundary condition that 
there is no net current into the nonactive wall (i =0) is 
satisfied by the electron contribution to current only. 
This is justified since it can be shown that the electron 
contribution is approximately three orders of magnitude 


greater than the ion contribution to the total current. 
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IV. RESULTS AND CONCLUSIONS 


A. GENERAL DISCUSSION 
Computer solutions are reported here for three cases 
presented in Table I. The active site voltage was 35 volts 


in each of the three cases. Included in this section are 


TABLE I: Computer Solution Conditions 


CASE No (wall) n; (wall) Ne (anode) an; (anode) 


i 10 0.0 J-matched 0.0 
II "float" 0.0 not J-matched 0.0 
ied "float" "float" not J-matched 0.0 


_— aE a ee eee eee ee ee ec ee ee a ee ee SO oo 


also some pertinent results obtained at 10 volts to illus- 
trate the effects the net production term has con the conver- 
gence to a boundary layer behavior for the species density 
and free stream electric field. 

The solutions are depicted as oblique presentations of 
the potential, species density, space charge, net produc- 
tion, electric field and Joule heating distributions over 
the "x-y" coordinates defined in the problem formulation. 
Also included are two-dinensional plots of various terms 


along a "cut" or line perpendicular to the electrode wall, 
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starting at the active site, with some starting at the wall 
along the inter-nodal line of symmetry (see Pig. 3). The 
solutions are organized as Case I,II and III results, will 
be presented as Pigs. 12,13 and 14 respectively at the end 
of this section. The discussion of these results proceeds 
With additional figures which are interspersed in the text 
appropciately. 

Table I presentS various electron and ion boundary con- 
ditions aiong the wall some of which are "floated." By this 
term is meant that upon the completion of each numerical 
iteration, the electron (and/or ion) density at the nonac- 
tive portion of the wall is re-established by assigning it 
the value of a near neighbor away from the wall. This pro- 


cess is illustrated in Fig. 4. 





DENSITY > 


x 
INITIAL CONDITIONS AFTER ITERATION "FLOATED" 


Pigure 4 


Description of "floating the Wall Boundary Condition. 
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The usual elliptic boundary value problem has well 
possed boundary conditions from a priori information. How- 
ever, the nature of the problem in this study includes very 
little information as to the correct species wall density. 
It seems reasonable to set these densities to zero along the 
inactive portion of the electrode as previously mentioned in 
Sec. II.B. In this study the wall boundary conditions are 
"“lossened" to facilitate "naturally occuring" and unpresun- 
ing species density profiles. The nature of the problem in 


effect, produces a self-consistent boundary condition. 


B. E@FECTS OF NET PRODUCTION 

Earlier work (Ref. 15] reported limited resuits for a 10 
volt active node site. Solutions for cases with and without 
the inclusion of the net production term have profound 
influences on the appearance of the charged species density 
profiles. Figures 5,6 and 7 illustrate the before and after 
effects of including the net production term. In both cases 
the wall boundary values for the species density was ini- 
tialized at the free-stream value (i.e., 1.0), then allowed 
to "float" as described previously. The addition of the net 
production term brings abont the boundary layer behavior of 
the densities of the electrons and ions is illustrated in 


Fig. 7. zfforts to attain converged solutions with the net 
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production term set to zero proved fruitless for active site 


potentials exceeding 10 volts. 





Ebectcon Density Ion Density 
(a) (D) 


Figure 5 


(a) Blectron and (5d) [on density perspectives with 
the net production term set to zero. (10 volts) 
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Pigure 6 


(a) =dhect ron and 


(b) Ion density eee with 
the net orctuction term couple v 


ae of olts) 


by 





SPECIES DENSITY 





O 
O I 2 i .4 ee .6 
DISTANCE (MM) 
@ = 10 volts ; 
i 2.35 10 ~ meter isothermal 
Fagus) 7 


Fluence of the net production term on the 
density’ profilas along 2 line perpendicular to the wall, 
ending from thé active node site, (10 volts). 
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oe eee USS@ON OF RESULTS 
1. Voltage Potential 

The voltage potential of the active node site ("the 
anode") in each of the cases presented was’35 volts. Fig- 
ures 12-a, 13-a and 14-a will clearly illustrate sharp vol- 
tage spikes at the anode, and the "Laplacian" nature of the 
field far from the anode. The potential along the inactive 
wall varied from 12 volts at the farthest position from the 
anode to 17 volts just next to the anode. Similar profiles 
of the electric potential along the wall for the other cases 
are presented in Fig. 8. Figures 12-1, 13-1 and 14-1 will 
show the voltage potential variation along a line normal to 
the wall extending from th2 anode site. 

Figure 9 illustrates a two-dimensional cut of the 
voltage potential along the inter-nodal symmetry line. Here 
it can be seen that the variation from the wall is shallow. 
It is important to point oat that the electron temperature 
is directly related to the gradient of the voltage potential 
through Eq. 5. This fact has certain consequences which are 


discussed in the next section. 
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2. Electron Energy 


The electron energy will be presented in series -b 
of Pigs. 12,13 and 14, as oblique perspectives. The anode 
site temperature for Cases I, II and III are 2.1, 2.1 and 
1.9 e¥ respectively. The description of the electron temp- 
erature is that of a linear function of the logarithm of the 
electric field. Any portion of the potential field that is 
"level" (grad Z =0) will necessarily generate a low or even 
negative (see Fig. 14-b) electron temperature. According to 
Fig. 2 and the discussion in Sec II.A.4, a negative tempera- 
ture would occur below an E/n of 0.04x107-29 Y-em2. This 
problem can be removed by setting a minimum value to the 
non-dimensional temperature field, nominally 0.13, corres- 
ponding to the minimum at the low end of E/n in Fig. 2. 

The parameters of the discharge (Sec. I.C) require 
an electric field of 2.7x105 V/m (or, logf{ BJ=5.43) at the 
free stream boundary, which corresponds to an electron tenp- 
erature of 0.905 eV (@=1.). The temperature fields in all 
cases presented do not quite meet the free stream boundary 
reguirements. The electric fields are insufficient to ren- 
der proper results to the electron temperature profiles. 
Thus the free stream temperature boundary requirements are 


"forced" by slightly inflating the influence of the electric 
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evelad. The effect of forcing the temperature field in this 
way can also be looked at as "sliding" the linear 
representation of the temperature curve to higher values at 
constant E/n, (see Fig. 2). Figure 10 illustrates a case 
where a low electric field, here 0.8 kV/cm, results ina low 
electron temperature projecting toward the free strean 
(9=0.62). The low temperature is adjusted to meet the 
desired free stream value (@=1.0) by multiplying the elec- 
tric field by a factor of 3.40. Pigure 10 also illustrates 
an alternate method of meeting the free stream conditions by 
truncating the variation of temperature with a minimuo 
value, namely the free stream boundary condition. 

Both of the "forcing" techniques are used, but for 
different reasons. As long as the adjustment factor does 
not exceed an order of magnitude change in the electric 
field, then it is not objectionable in light of the simpli- 
fied representation of the electron temperature. The 
adjustment requirements of the three cases reported in this 
study are between 3.0 and 5.0, with Case I using a factor of 
10 to see the effects on the Joule heating. In Case [I the 
non-dimensional temperatur2 is artificially high at 1.4 eV 


at the free strean. 
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The electric field near the inter-nodal symmetry points 
of the wall are characteristically 2 orders of magnitude 
less than that of the free stream conditions. The resulting 
very low, even negative temperatures in Case III, Fig. 14-b, 
is "patched" using the truncation method wae ogi dis- 
cussed. Whether these low temperatures and the associated 
steep derivatives are physically reasonable, merits further 
investigation. 

The temperature fields generated and presented as 
Figs. 12, 13, 14 series -b are the final result depictions. 
However, the temperature fields used in the solution process 
were adjusted and truncated as previously discussed to main- 
tain solution convergence. An example of a "corrected 
temperature field is presented in Pig. 14-c. 

3. Electric Field And Joule Heating 

The electric field and Joule heating distributions 
are generated from the primary voltage and electron density 
solutions, much in the manner as the temperature distribu- 
tion was determined. In fact, the Electric field and the 
Joule heating distributions are similar in appearance to the 
temperature distributions. These distributions play no 
direct part in the soneion scheme and thus being ancillary 


in nature were not "forced" to meet any a priori conditions. 
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The Joule heating term can be written as; 
JoeB = [DML KT / Ay? Jel VG -Vne 1°74 

and thus is a rather straight-forward function of @, ng, 
and n;. The ion contribution is shown to be negligible, 
Since D,>> D;. The Joule heat at the anode site is of order 
108 and is 100 times that of the free stream values. Joule 
heating effects are negligible in the overall problem, as 
discussed in Sec. II.A.3, obviating the reguirement for an 
energy equation. 


4. Species Density Distribution 


== 2 ae ee eS 





The solutions for the charged species density dis- 
tributions are presented as oblique perspectives in Figs. 
12,13 and 14, series -e (electron density); series -f£ {ion 
density); and series -g (space charge density). Addition- 
ally the variation of these distributions normal to the wali 
along the active node symmetry line are presented in the 
series -j and -k figures. The oblique perspectives of the 
densities are reversed for clarity. 

The electron and ion densities at the nonactive wall 
are a subject of interest in the literature, to which little 
information has been gathered on a firm basis. The results 
of this study suggest a slightly positive space charge at 


the inactive wall. Scrutiny of Case III results reveais, a 
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2.2x10!© m-3 positive space charge number density, (or 3.5 
milli-coulombs). Recall that case III involved "floating" 
both of the charged species densities at the wall. In fact, 
close inspection of Case I and II reveals a positive space 
of similar magnitude but displaced slightly from the wall. 
The wall in these cases was set with a negative space charge 
as the boundary condition. The positive space charge along 
the wall is most predominant in the vicinity of the inter- 
nodal symmetry line, and whether the low temperature values 
in this region have significance, remains to be determined. 

The active node site (anode) electron density in 
Case I (2.4x10!7 m-~3) was determined so as to match the free 
stream current to the current constricting at the anode. [In 
the results of Cases II and III the anode electron densities 
are chosen so as to provide a smooth variation extending 
outward from the active site. 


5S. Tonization/Recombination 





The net ionization term is utilized in all cases; 
and was found necessary for (1) establishment of the charged 
particle density distributions, and (2) solution convergence 
for sheath lengths spanning more than one computational node 
spacing. The net production term is depicted in Figs. 12,13 


and 14 series -h. There is net ionization taroughout the 
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solution domain. However as Case III demonstrates, there is 
a factor of 5 orders of magnitude greater ionization within 
the sheath region surrounding the anode site than at the 
free-stream edge. Cases I and II show ionization along the 
nonactive wall, simply because the ion density is zero, eli- 
Minating Significant recombination, (see Eq. 4%). 


6. Current Description 





Figure 11 is typical of the constricting nature of 
the current stream lines. The potential, electron density 
and the electron temperature are used to determine the "x-y" 
components of the current density. Using a "nearest neigh- 
bor" influence principal, a local derivative can be deter- 
mined; from which a streamline can be "marched-out"™ in small 
increments. The profiles of all three cases are similar, 
the fact being that the current matching discussed previ- 
ously will not affect the current streamline profiles. Fig- 
ure 11 clearly demonstrates the satisfaction of the current 
constriction requirement necessary to permit solutions to 
the system of equations. 

The current into the anode site in Case I is matched 
to the current entering at the free-stream boundary. This 
specifies the electron density at the anode site. However, 


in Cases II and III the current at the anode was only three 
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times the free-stream current. The current matching problen 
does not deserve serious complaint, since the problem 
description is simplified tc two-dimensions and the current 


constriction is only required to satisfy continuity. 
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Two-dimensional depiction of the current stream line 
ee on ecard the active node site. 
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Figure 12-b 
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Oblique perspective of 
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Pigure 13-a Figure 13-b 
Obligu2 perspective of Iblique perspective of 
Pocgencial Field (CASE ITI) Temperature field (CASE ITI) 
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V. CONCLUSIONS AND RECOMMENDATIONS 


= es 42 SP 2S es EP a =e EO SE ce TS a SS LS ee eee es oe 


A> REVIEW OF ASSUAPTIONS 

In the development of this work several assumptions were 
made which the solution results later verified. The assump- 
tion of negligible convection effects is straight forward. 
The charged particle density boundary layer evolution time 
turns out to be small in relation to the pulse period; thus 
enabling consideration of the problem as a steady-state aft- 
er-glow probiem. The fact that the Joule heating of the gas 
is shown negligible, frees the formulation from an elaborate 
energy description. Since the electron temperature is a 
Simple function of E/n, the prodlem can be further simpli- 
fied. The solution proceeds with nominal values of the par- 
ameters of the discharge as discussed in Sec. II.B. 
Finally, the importance placed on the net production term is 
verified by the profound changes in the results including 
the fact that "ng" is necessary for the establishment of the 
charged-species boundary-layer behaviour. The analysis of 
App. 3 shows that an assumption of "frozen flow" would in 


fact pe preaature. 


78 





BB. PHYSICAL CONCLUSIONS 

The present model of the sheath and ambipolar region 
evolves from the assumptions of steady state, uniform ion 
temperature and constant diffusion coefficients. A summary 


of important results is presented in Table II. 


TABLE II: SUMMARY OF IMPORTANT RESULTS 


SUMMARIZED CASE CASE CASE 
RESULTS I {I Ieee (units) 
Sheath length 35 3.5 37 10-5 o 
ambipolar length 4.5 55 Sie 10-5 no 
Eo (free-strean) Teo S09 6233 10* V/n 
= (2neode ) io Des Zee 106 V/n 


The results show that a sheath develops as anticipated 
from theory and the proposed assumptions. The sheath char- 
a@eteristic length is consistent with that predicted by Debye 
shielding. Debye shielding of a flat plate predicts that 
the sheath length "A," is 4.4x10-5 meters. Por a point (or 
spherical) shielding the sheath length is 3.1x10-5 meters. 
The present study involves the shielding of a point voltage 
source situated on a surface charged plate. Since both two- 
and three-dimensional effects are inherent in the formula- 


tion, the resulting sheath lengths are quite reasonable. 
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AS Table II shows, the size of the ambipolar region is of 
the same order as the sheath size. Here the ambipolar 
length is defined as the length extending from the edge of 
the sheath to the point at which 99% of the free-stream den- 
Sity is recovered. FPinally, the sheath can be self-generat- 
ing from a consistent set of equations without an overall 
energy equation. 

The resulting electric fields at the free-stream boun- 
dary do not guite attain the desired 2.7x105 V/m presented 
in the problem formulation. These results fall short by the 
same factors required for the adjustment of the electron 
temperature field to meet the stated free-strean boundary 
reguirements. It must be pointed out that the electric 
field is the reguired boundary condition. However, this 
boundary condition could be made to yield no reasonable 
results, and so the chosen boundary condition is © = 0. 

The resulting electric fields at the active site in each 
of Cases I, II, or III was of two orders of magnitude higher 
than that of the free-stream. These electric fields are 
below that needed for breakdown of Nitrogen at one atmo- 


sphere, which is determined as, [R2f. 18]; 


Zh = (E/p)Pp = (350 V/em-torr)x(760 torr) = 2.7x107 V/n 
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The requirement of strictly matched currents produces 
steep derivatives and inflections of the density profiles in 
the vicinity of the anode. A noticeable "spike" is evident 
in Fig. 12~j for example. Case II and III anode electron 
densities were established on more aesthetic grounds, where 
snoother appearing density profiles were generated. In 
these cases the anode current exceeded the free stream cur- 
rent by only a factor of 3. 

The nonactive regions of the electrode (wall) far from 
the active site reveal a small positive space charge den- 
sity amounting to approximately 2% of the free-stream den- 
sity. The negative space charged sheath is 46% to 56% of 
the free-stream density. 

Previous work [Ref. 12] of a Similar problem reported no 
reasonable results for the case of a catalytic wall 
(n, =n; =0). The results presented in this study demonstrate 
wall conditions of n-= 0. and ne < 1. Here, a noticeable 
positive space charge coexists with the large negative space 
charged sheath, without large inflections in either the vol- 
tage potential or densities. 

The Case III study has both electron and ion wall condi- 
tions floated. The shape of these density distributions is 


quite symmetrical about the anode site in all directions. 
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One needs little imagination to visualize that the density 
profiles can be "wrapped" around the anode completely, in 

effect treating the anode as a whisker or protrusion into 

the plasma. However, the imagination must be tempered by 

the fact that it can only be a two-dimensional protrusion 

(e.9., a2 *blade). 

The effect of the net production term has been previ- 
ously discussed as one of the most important results of this 
study. The term is necessary for inclusion in the problen 
in order (1) to attain wall electron density values of less 
than the free-stream, (2) to produce the boundary-layer 
behavior of the species densities,and (3) to obtain con- 
verged solutions with voltages greater than 10 volts or con- 
putational node spacings of any order less than that of the 


sheath length. 


C. NUMERICAL CONCLUSIONS 

The Modified Newton-Raphson (MNR) method, presented to 
solve the particular set of coupled non-linear partial dif- 
ferential equations in this study, 1s a general method that 
can be applied to more sophisticated systems. The defi- 
ciency of the quasi-Jacobian method (discussed in Sec. III.C 
and App. C.) in dealing with non-linear terms is circun- 


vented by an alternative linearizing scheme. Convergence is 
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rapid; however, the starting solution must not exceed 10 
volts, must have zero net production, and must have a cons- 
tant electron temperature during the first 2-3 iterations. 
After the solutions take their initial form at ten volts, 
the potential field is "amplified" and reinstated as a new 
Starting solution. The final stages of the converging solu- 
tion are carried out under the full influences of (1) float- 
ing wall densities, (2) variable electron temperature, (3) 
net production and (4%) an active nod2 site electron density 
being solved for through current matching. 

The profiles of the space charge and net production dis- 
tributions demonstrate erratic behavior at the first few 
computational nodes from the active site. This is believed 
to be caused by the linearization of the governing non-li- 
near equations and has the appearance of truncation error 
when observed on an iteration to iteration basis. However, 


the remaining regions of the profiles are smooth. 


D. RECOMMENDATIONS 

This work succeeded in solving the problem in a two~di- 
mensional environment. This fact necessarily limits the 
quantitative results since the current density decreases in 
only two dimensions. With the computational power of recent 


coaputers, it should prove worthwhile to make the next 
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advance as a three-dimensional problem rather than elaborate 
in the two-dimensional case. A three-dimensional scheme 
would render more realistic results of a rere nature 
because active spots are more likely to exist than active 
lines and because 3-D geometries require less voltage for a 


gzven current. 


— 


Specific recommendations other than the above include; 


e Account for the ion contribution to the current 
boundary condition into the wall. 


e Investigate the physical See Cohn e of the low 
electron temperature near the wall far from the 
anode site. 


e Investigate the physical implications of the elec- 
“ron density "spikes" prevalent near the anode site. 


° ra the wall temperature a valid boundary condi- 
ion. 


e Use a free-strean ote condition of . 
Eoo =2.7x105 Vym, rather than zero potential. This 
would introduce desired results for pumping lasers. 


e Include a more sophisticated electron temperature 
model, or verity the present one. 


e Increase the computational mesh density. 


e Study the general application limits of the MNR 
method. 


e Investigate "n " in the role of making one-dimen- 
Slonal Solutions possible (1.e., alter the reacting 
rates)... 


e Investigate whether the low electron temperatures 
and assOciated steep derivatives located along the 
wall between active sites have any fee implica- 
On or result merely from inadequate problem for- 
nulation. 


e Consider the electron-beam ionization in the solu- 
tion to produce more realistic effects. 


84 








APPENDIX A 
CHARACTERISTIC PARAMETERS 
A. SHEATH LENGTH 
The extent of the sheath is one of the most important 
G@arcteristic lengths in this study. It is within the 
Sheath that most of the potential drop occurs for a short 
discharge. FPortunately, the sheath length can be estimated 


rather easily from Poisson's equation, (Eq. Ai) 





v'¢ = €- Me-ni) (44) 
Let; 
A - A 
P= ate! "i a 5 Ns = Ms Neoo (A2) 


#here ")," is the characteristic sheath length; the sub- 
sctipt "o" indicates the value at the electrode and " " the 
value at the undisturbed plasma. After gathering the dimen- 
sional quantities on one side, Eq. A2 becomes; 


ee 
= e New As -|\Ne-Ni (A3) 
V4 151 | e \ 


The characteristic sheath length " s'' defined in this 
fashion, is used to non-dimensionalize the "del" operator in 


all directions. That is, independent of the dimensionality 
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of the problem, the sheath is effective over the character- 
istic length just defined. Now an estimate of the size of 
the sheath can made by setting the grouping of dimensional 
quantities in Eq. A3 to be of the order one, namely; 


A= isa] = Ap sig" (A4) 


C Ness K C90 


where Ap," is the Debye length in the undisturbed plasma. 

Anticipated values are now established for the sheath 
length. Take for example a 35 volt anode potential, anda 
Charge species density of 10147 to 10!9 n-3. These parame- 
ters define a sheath length of 1.4x10-* to 1.4x10-S meters. 

AS can be seen, a ten-fold increase in voltage would 
change he by a factor of 3, rendering the nondimensionali- 
zation of the problem relatively insensitive to the chosen 
value of the voltage potential. That is to say the value of 
the sheath characteristic length as defined may be consid- 
ered a resonable estimate for the charge species densities 
examplified. 

The important implication here is that extreme care must 
be devoted to electrode surface preparation to preclude 
whiskers or flaws of a size comparable to nA, These pro- 
trusions would produce a focusing of the electric potential 


to create extremely high electric fields, which become sites 
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for initiation of breakdown. Protrusion growth has been 
observed to form in a time on the order of milliseconds. 
Many conditions are favorable to protrusion growth, depen- 
dent on the pre-breakdown and breakdown current profile, 
electrode composition (especially in the case of a ther- 
mlonic cathode), surface purity and surface preparation 
(Ret. 21]. Some possibilities are; 


e Electrode plastic flow or melting in an electric 
field due to micro-particle bombardment. 


* Stickmmag of Wicroparticlées. 


e Electron beam induced protrusion growth. 


Be. DIFFUSION LENGTH 


If diffusion is a significant factor in this study then 


its contribution to the current would be; 


a ONe _ on. egne 
J, = Dees D: a 


or, since Dy > D-, and using Einstein's relation, or 
O/ he =kT, {23 


J x 204 5 = = pe kT, St 


making the equation unit free (through fractional analysis) ; 





JL an CNG opper \t\ (A5) 


@ De Ney o x 
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If the right side of Eq. AS is of the order one, (the nondi- 
mensional density gradient is a maximum of 1) then "L" may 


be estimated as; 


Ls eDeng /J = 3.7x10-¢ on 
for the values; 
ne 10'S a-3 
De 0.0665 m2/(Vesec) 


Je =Heeng VF = Hercongen,e(E/n) = 2.88x103 A-n-3 
Since &> L, then any role attributed to diffusion 


shculd be visible within the sheath. The E/n within the 
thin sheath is of course greater than the E/n at the free 
Stream (1029 Vem2), which decreases the local value for "L" 
ead further. 

Another approach is tc assume a one-dimensional sclution 
to the charged species conservation equations, which can 
only be done outside of the sheath, i.e., in the ambipfolar 
region. The ion conservation equation may be written; 

n = Veg = eVe(YnVP - Dvn] (A6) 
Since it can be shown (next section) that the production 
term can be approximated as; 
nxn -Xn2 
then for the stationary condition of no net production; 
Y= OLN 


Nondimensionalizing fq. A6, it can be shown that the 
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Characteristic diffusion length "L" is; 
L = YDACN = 2.x10-© meters 
which is of the same order of magnitude as the characteris- 


tic sheath length. 


C. ILONIZATION/RECOMBINAITION LENGTH 

Analysis of the production term (Zq. 4) reveals that for 
Nitrogen the three-body (electron) collisional recombination 
cate (X%,) is small in relation to to the other terms, (see 
App B.-). In any case the net production term is; 

Ne = no- Cnn, ~ %n en, Dy (AT) 

In the steady state the net production term is zero, and 
as is done in Appendix B, the relation becomes; 

De = (&,+ Ong) o(n%- ng) ong 
Substituting the value of 0.17x107-!!1 m3/sec, and (nt - No) & 
1018, we can eStimate the characteristic production tine to 
be Tp = 5.88x10-7 sec. 

To establish a meaningful characteristic length, the 
average thermal velocity will be used since it represents 
the principal transport velocity of the mediun. 

fp=Tpeve = 5-88x10-703.2x102 = 1,88x10-% a 
The characteristic length for net production is a factor 


of four larger than the sheath length ee 4.4x10-S5 m). The 
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conclusion which can be aade from this (fractional) analysis 
is that the net production is significant within the overall 


region of interest, including the sheath itself. 


D. FROZEN FLOW 

When a process relaxation tine ae is long compared 
with the time that a fluid particle spends in a given 
region, the flow is said to be "frozen" in this region. By 
this term one means that a reaction process may be 
neglected, as far as influencing the flow is concerned, and 
the species involved treated as chemically inert, [Ref. 22}. 
In the opposite limit, where "%" is very small relative to 
wk the characteristic flow time, the flow is said to be in 
feguilibrium." for the recombination rates given in this 
section, wpe = 9 -Xx10-7 and ne = 4.0x10-?7 seconds. Here 
flow time is defined as the sheath Characteristic length 
G@ivided by the drift velocty. Thus the flow is not “frozen" 
as far as chemical reaction, but we must consider that the 


electron and ion species are not in thermal equlibriun. 
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APPENDIX B 
DEVELOPMENT OF THE PRODUCTION TERM 

The net production term is defined as the difference 
between lonization and recombination rates. Whenever a neu- 
tral atom or molecule loses one or more of its electrons, 
ionization has ocurred. Recombination is the process wher- 
ein positive ions and electrons combine to form neutral 
atoms or molecules. It 1S simply reverse lonization. If 
the net production is zero then ionization and recombination 
are equal. The net production term is the right side of 
Eqs. 2 and 3, the expression for which is deveioped here. 


(1) Dissociative - the electron recombines with a molecu- 
lar ion (AB)+, and the recombination energy goes 
into dissociating the molecule and iacreasing the 
kinetic energy of the resulting products. The pro- 
cess is described as; 

e + (AB)+ — » A + B 

(2) Three-body (electron) - the energy released in the 

recombination process is carried off by a third body 


wnich is an electron. The process is described as; 


a+ At + So ——%» A t+ 
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(3) Three-body (heavy) - the recombination energy is 
carried off by a third body which is a heavy parti- 
cle such as an atom or molecule. The process is 
described as; 

e+At +B —»> A +B 

This study does not consider other recombination mechan- 
isms such as radiative or dielectric recombination, since 
they are generally of low significance in relation to the 
just described processes. The ions are limited to singly 
positive ions. No account is made for intermediate or 
excited states as this reguires a detailed balancing which 
is beyond the scope of this work. Additionally, no account 
is made of attachment. 

Rewriting Eq. 4 for simplicity of discussion; 

De = Y- De - Ono n; - HN *n; — Xe ne in (B1) 
where each "gw." refers to the recombination coefficients for 
dissociation, three-body (electron) and three-body (heavy), 
respectively. 

The recombination coefficients depend upon the ionic 

species involved. Some general but representative values 


(Ref. 22] follow; 


”® = 10-13 m3/sec 
Wy, = 9x10-'* n3/sec 
Xx = 2x10-13 m3/sec (air) 
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Hinnov and Hirschberg [ Ref. 23] developed a relation for 
three-body (electron) recombination coefficients as a func- 
tion of electron temperature (Te) and density, vaiid for 
energies less than 0.25 eV. 

x, = ee xa0—2°ne Test 7< y37sec 
where the electron temperature is °K. Gurevich and Pitaen- 
Ski [Ref. 24] extended this development to obtain; 
My= 2.13x10-29n Te-9/SlnA 03/sec 
where; 
InéA = 1.24x107 (Tp 3/n,) - § 
At 1.0 eV these coefficients are 5.6x1072! and 1.05x10-19 
m3/sec respectively. A simplification of the later expres- 
sion for @ is valid from energies of 0.1 to 10.0 eV and 
a so-* © a; 
AX =0.032xT-*4- 3 

poewn (Ref. 18: p.197] gives a value for the dissocia- 
tive recombination coefficient; 

&, = 1.5x10-12 mn3/s 

The ionization coefficient can be determined from con- 
sidering that there is no net production in the free strean, 
i.e., far from the electrode. Her2; 


0 =4N -XKN2 - XN3 - 4N2n, (B2) 
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where "N" is the stationary state. Solving for and sub- 


Srieimeng into Eq 31 gives; 


ne= Np (N-n; ) (%j+%No) +O,neTe9/20 (N2-n, n;) (B3) 
where; 
@, = 1.5x107-t2 (from Ref. 18) 
, De= 1.09x10-290 (from Ref. 23) 
&ng= 2.0x10-13 (from Ref. 22) 


The net production term is the relation, suitably non- 


dimensionalized, that appears in the computer subprogram 


mao DDE. 
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APPENDIX C 


NUMERICAL SOLUTION METHOD 


A. GENERAL DISCUSSION 

There are two basic approaches used in solving a nonli- 
near syst2m of equations, such as possed in this study. A 
given method can usually be classified as one of the ascent 
(or descent) or more commonly refered to as gradient meth- 
ods, or as one of the Newton methods. The method applied in 
this study and by Dolson [Ref. 12] 1s one of the Newton 
methods. A brief review of both methods will be presented 
to facilitate development of the particular method employed 
in this study. 

Let the nonlinear system of equations be written as; 

£; (X,eXqeee een) = 0 


ie (C1) 


| 
>) 


Ey (X, oXgaee eX) 
or in vector form for convenience; 
See “~~ 
£(x) = 0 (Cz) 
In the descent method the problem is possed in such a 
way so as to form a function which is zero at any selution 


and positive otherwise. A typical such function would be; 


gS 








n 


iGo ces xr) =il (C3) 


so that the problem is reduced to finding a minimum of @, 
ahich also happens to be the value zero. In solving prob- 
leas by this method, an initial guess to be solution must 
be provided, call it x€!)2, Then a direction d is determined 
sO aS tO point t9 a better solution, x2). This direction 
is often chosen to be the gradient of (x). The operation 
-grad @(x) is the direction in which @(x) decreases most, 
and is locally the best direction in which to reduce @(x). 
d = -grad (x) 

so that; 

x€2) = x€1) - aegrad g(x¢?)) 
where “a 1s chosen to be an appropriate step size. A prob- 
lem is encountered at local minima of @(X), which needs to 
be avoided. One must also use the step "a" in such a manner 
as to control the acceleration along slowly or rapidly vary- 
ing @(x). It may be that by the nature of (xX), the gra- 
dient may have to be approximated by finite differences, an 


approach which offers a flexible recourse even when not 
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The Newton metnod involves approximating the nonlinear 
problem by a suitable linear problem involving derivatives 
of the nonlinear equations evaluated at x€!). The linear- 
ized problem is simply a set of linear algebraic equations, 
the solution to which furnishes an increment that when added 
to the solution x!) will yield x2). This method is a gen- 
eralization to n-dimensions of Newton's method for evaluat- 
ing the zeros of a function. 

Each function fj (x) is approximated by the initial two 
teras of a Taylor series expansion about the solution vector 
x1); 

Gey = ty eh) + | Jeo pet x = x6) | (C4) 

where "j]" is the familiar Jacobian matrix, 

n 

j= 3f4(x'” ) (C5) 

med Xj 

j=l 
The right-hand-side of equation c4 is the linear vector 
function (tangent hyperplane) of x which best approximates 
the nonlinear function f evaluated at x€!9, We set this 
right-hand-side equal to zero and solve for x, which is an 
improvement from x6!)2, We call this improved solution vec~ 
tor x2), and from equation C4 we Jet; 


x€2) = xC1) = J(eO 19) m2 ef (x) (CS) 
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The Jacobian need not be known very accurately, and often is 
left unchanged for several iterations since its evaluation 
often entails penalizing computer time. Thus the method 
(slightly modified as described briefly in section IIIZ.B) 
used by Dolson (Ref. 12], is referred to in many texts as 
the Newton method or as the Newton-Raphson method [{Refs. 25, 
26 j. 

The Jacobian matrix of coefficients for elliptic prob- 
leas is typically a sparse banded matrix, to which there are 


Many computer subroutines available for solution. 


B. MODIFIED NEWTON-RAPHSON METHOD 

The particular solution technique employed in this study 
merits a detailed discussion, since it is not common in the 
literature. Bailey and Touryan {Ref. 27] use a similar 
method and provide some description of its development. The 
method is Similar to other Newton methods however, the Jaco- 
bian matrix is replaced by a matrix of coefficients devel- 
oped from a sequence of linear differential equations which 
approximate the system of nonlinear equations. This method 
Will be referred to as a nodified Newton-Raphson method or 
MNR method for short, and will be presented in general using 


a specific example. 
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Take Poisson's equation given eariier and repeated here 
nor siaplicity; 
V?$= sc, (Me -n,) (C7) 
The known Solution variables are; 
PCtr, Nett), and n; ¢1) 
In general tne £Eirst step is linearization of the equa- 


jem, in) thas case Eq..c7? linear, but it will serve as 


la 


an example. We assume that the suceeding solution variables 
at iteration-(2), namely z°2), can be approximated by the 
solution z61!) plus some small change Az‘ 1), 

QAzc2) = 2741) + z(1) (C8) 


where; 
cK) = CK) CK) Ck) 
< = (¢ re , Ql; ) 


It remains to determine 42°12. Substitution of this 


' Cs] q@) 
ve” 2 oe + vag’ = Gig n°) + (Ne - AN, ) 


rearranged; 


' < y 
vagt?c,(ane'- Ani”) eee [we ve Z, (ne-n3 )| (C9) 


She, 








The right hand side of Eq. C9 is a function of the known 
Solution vector z¢1), We are applying this relationship at 
each computational node. A second order partial differen- 
tial equation can be approximated by a second order differ- 
ence equation or "3-point formula," so Eq. C9 will become a 
system of linear algebraic relationships involving near 
neighboring values of the solution vector z‘!). The left 
side of Eq. C9 is composed of a system of coefficients 
expressible in terms of z(1!) times the unknown (perturba- 


tion) vector 42z!). So we now have; 
A(z61)jeAze1) = -P(z1)) (C10) 


Here the matrix of coefficients "A" is a sparse matrix 
With a relatively narrow band of non-zero elements about the 
Main diagonal, provided that the "natural" order of the 
resulting linearized equations is preserved. There are 
standard computer routines available which take advantage of 
this banded structure, such as LEQT1B, an IMSL routine. 


It follows that since Az(1) = z¢2) - 721); 


zo2) = 2O1) = FA(zO1)) J-teoPr(z1)) (C11) 
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For the case of a nonlinear term such as: 
No oF 
° 3x 
there is no additional difficulty. When linearized this 
term becomes; 


Ne2F sAnedd JAD | Ane JAG 
ed re tee ~ e x (C12) 


We neglect high order terms created by the inner products of 
the elements of Q2!), such as the last term in relation 
C12. The first term would be a part of the function F(zC!?) 
in the right side of Equation C10 , and the remaining two 


terms would be included as part of the expression: 


A(zCtdype Azer), 
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APPENDIX D 


DESCRIPTION OF COMPUTER PROGRAMS 


This appendix serves as a user oriented guide to the 
computer programs listed in App. E. The programs are listed 
in the order used to generate solutions to the problem for- 
mulated in this report. The output consists of printed 
results for each iteration and punched output to provide a 
permanent file copy, which is to be used as input to the 
three plotting programs also listed in App. E. These plot- 
ting programs use the Versatec plotter and associated soft- 
ware, and represent the only facility not commonly available 
elsewhere. The main program SHTH as presented, is config- 
ured to also produce output at a monitor terminal as used on 
CP/CMS (IBM 360) or MVS (IBM 370). The output in summary is 


Via three devices; 


PaNLCe & = typically a CRT 
DEVICE 6 - typically a printer ; 
DEVICE 7 - typically a bunch or disk output 


These can be manipulated in convenient configurations to 
give a full or partial (sampled) display of the solution at 
each iteration, or merely to display convergence information 


in 4 Simple format. = 
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The SHTH program is the main program which together with 


its 9 subroutines and 1 IMSL subroutine solves the non-di- 


mensionalized forms of Eqs. 1 through 4. The input to the 


program 


] 
3 
1 


consists of control cards and data cards, namely; 


control card | 
or 104 initial solution data cards 
control card 


The control card is described here in terms of the Fortran 


Variable 


ITAIN 


ITNAZ 


ITNDT 


ITEAP 


ITBC 


LADJ 


LPOS 


PRT 


NAME 


Names thereon. 


- is the iteration counter (ITER) starting value. 
If ITHIN = 1, then the program will read in the 
next 3 cards as data which’ contain the short forn 
Starter solutions. If ITMIN is greater than 1 

then the program will read in the next 104 cards as 
data, cepresenting a detailed starting, solution. 
This detailed starting solution is typically the 
Output 29f a previous run. of the program. I£ ITHIN 
= 00, then the program will not read in any data, 
being thus informed that the problem is over. 


- is the last iteration value. The program ceases 
1terations at this point and outputs the solutions 
and auxillary information. 


- enables the employment of the net Peo ae term 
(RTSIDE) when ITER is greater than ITNDT, otherwise 
the net production is set to Zero. 


- for ITFER less than ITEMP, the electron tempera- 
ture is isothermal, otherwise it is a function of 
E/n. (ses SHTH statements 899-999). 


- for ITER greater than 5 and less than IfTBC, the 
Peery conditions for the Species 15S 
oated. 


- 1 electron temperature is calculated 
- 0 electron temperature is isothermal. 


- 1 output is punched on device 7. - 9 output is 
not punched. 


- sets the mesh spacing equal to H*ZAMS (sheath 
length). 


-~ is the adjustment factor to the compere lube equa- 
tion (see Séc. IV) and is used when LADJ =1 and |. 
ITER is less than ITEMP. The adjustment factor is 
not used when LADJ is_set to zero. When the temp- 
2rature is coupled and LADJ = 1, then the program 
calculates its own adjustuaent factor. 


- is reserved for user labelling and identification 
which is used at output. 
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The solution set is; 

PH voltage potential {volts). 

NE electron density (10!8 7-3) 

NI ion density (1018 m-3) 

TH nondimensional electron temperature 

The short form starter solution is illustrated in App. 

E. It consists of PH, NE and NI values at the anode, wall 
and equilibrium boundaries. The long form starting solution 
consists of PH, NE, NI, and TH distributions, each one 
represented by 25 data cards and a label card. Each card 
has 7 values. (see format-803). 

typical iteration series would start out for about ten 
izerations from the short form data with ten volts at the 
anode and other appropriate boundary conditions. The ini- 
tial soluticn converges very well for # = 0.5, the net pro- 
duction term coupled after 2 or 3 iterations, and the wall 
conditions not floated. This output is then used as inout. 
The voltage is amplified to say 25 or 35 volts and a new 
control card is structured to couple the temperature or wha- 
tever else is desired. The H value may also be changed to 
reflect a smaller mesh spacing. It must be kept in mind 
that too many changes in One solution run can be catas- 


trophic. 
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Bach iteration of the solution by the SHTH program 
accomplishes the following; 

(1) floats the NE, NI wall values 

(2) sweeps through the A-matrix and C-vector setting 
up the equations that satisfy the boundary condi- 
tions and the internal mesh conditions 

(3) couples the electron temperature equation 

(4) couples the net production tern 

(5) calculates the anode electron density required to 
Satisfy current matching btween the anode and the 
free-stream. 

The system of equations has the form presented in Eq. 
C10 in App. C, and is solved by use of the IMSL routine 
LEQTIB. The solution vector is then updated from the return 
information contained in vector C. 

There are two solution convergence nonitors. The first 
is the sum of the squares of the differences (spectral norm) 
or the -F(z‘1)) values from Eq. C10. The other is the spec~ 
tral norm of the Az(t) upon return from LEQT1IB. Both of 
these are available at the video terminal (device 4) and 
printer (device 6). 

After solution update, the elctron temperatur2 is calcu- 


lated and adjusted per the control card options. [It should 
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be pointed out that on the last iteration, the temperature 
distribution is automatically calculated for output, but is 
not used in the scheme, unless otherwise coupled per the 
control card. 

The matching of the anode current to the free-strean 
current is accomplished by solving for the anode electron 
density using subroutine ANODE. The call to this subprogran 
will require a change to the SHTH program as it is presently 
Overridden with a comment card, (prior to statement 9000). 

The final outputs are straight-forward. Subroutines 
JWALL and JOULE are called in this phase of the program, in 
order to compute the current into the wall and the Joule 
heat distribution. The electric field, net production and 
space charge density distributions are also calculated and 
provided as punched and printed output. 

The punched output is used as input for the STREAM, PICT 
and PLOT programs. These programs use the Versatec plotter 
to provide a presentation of the current stream lines, an 
oblique perspectives of the solutions and plots of the solu- 
tions as they vary normal to the wall extending from the 
anode site. The details of these programs are combersome, 


and it is wise to leave them to the user for inspection. 
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The first two distributions in the punched output are 
used in the STREAM program, requiring that the nonpertinent - 
data be removed. The next 8 distributions are used for the 
PLOT and PICT programs. The PICT program has an option card 
filled with ones or zeros in the first eight positions. 

This option card enables the user to elect which distribn- 
tions to plot. The plots provided in this report as Figs. 
12 through 14, were produced using the PICT and PLOT pro- 


grams. The programs are simple and speak for themselves. 
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APPENDIX E 
COMPUTER PROGRAM LISTINGS 
This appendix contains a listing of the programs used in 
the solution of the equations presented in this work 
(namely: Eqs. 9,10 and 11), and the auxiliary programs that 


provided the plots for the report. 


SHTH (MAIN PROGRAM) -~-------------------- 109 
ee ne Fee oe ee COs se fo 
gt ae eo ee ee a. eee eee oe 45g 
mere — Seat ne ae 
Pane = Soe ee ee ee 
RR 2 a ee ee ee 133 
Wee eee Gea ee ea 
Wet = Sa an oe eee 136 
ae See -. Se See ee 137 
OUTPUT ------------ OES ee ei, ee eas er 138 
PLOT --- ~--- -- -- - = -- - - -- - - - = - - - = - - - - - --- = 439 
ee a eee ee eon a tes naa 
Ries Ses ee See kook k em 
INITIAL INPUT (sample) ------3~-3----------- 150 
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